function GSSOut = gs_itera(GSS,d0)
%gs_itera - Iterations for strain computation
%
%       GSSOut = gs_itera(GSS,d0)
%
% This function computes the strain rate field on the basis of grid nodes
% and velocity vectors in the struct variable GSS and add to GSS, leading 
% to GSSOut, these fields: 
%
%   MaxStrainMatrix     : maximum strain matrix;
%   MinStrainMatrix     : minimum strain matrix;
%   AzimuthStrainMatrix : matrix of azimuth, i.e. angle between north 
%                         direction and maximum strain vector;
%   MaxStrainMatrixError: error on maximum strain;
%   MinStrainMatrixError: error on minimum strain;
%   NumberMatrix        : numb(m,l) is the number of experimental points 
%                         whose distance from the grid node (X(m,l),Y(m,l) 
%                         is smaller than or equal to d0;
%   ComputationFlag     : is 1 in the grid nodes where the strain is 
%                         actually computed, 0 otherwise;
%   SignificanceFlag    : flag on significance. It is 2 for a grid node if 
%                         all the three sectors with 120� angular width 
%                         120� around this grid node contain at least one 
%                         experimental point whose distance by the node is 
%                         lower than d0 (HIGH SIGNIFICANCE GRID POINT CASE).
%                         It is 1 if two sectors contain at least one point 
%                         as above (MID SIGNIFICANCE GRID POINT CASE).
%                         It is 0 othewise (LOW SIGNIFICANCE GRID POINT 
%                         CASE). Obvioulsly, FLS is 0 if FLC = 0;
%   UTranslation        : matrix of velocity/translation East components;
%   VTranslation        : matrix of velocity/translation North components;
%   UTranslationError   : error on UTranslation;
%   VTranslationError   : error on VTranslation;
%   Rotation            : for each grid node, it is the non-zero component 
%                         of the rotation antisymmetric matrix);
%   RotationError       : error on rotaton.

% G. Teza, 2006, 2022.

Options = GridStrainOptions; 

magnitf = 10000;

[X,Y] = gs_gridExtra(GSS);

format long e

szx = size(X);

EMAX  = zeros(szx); 
EMIN  = zeros(szx); 
PHI   = zeros(szx);
eEMAX = zeros(szx); 
eEMIN = zeros(szx);
numb  = zeros(szx);
FLC   = zeros(szx); 
FLS   = zeros(szx);
U0    = zeros(szx); 
V0    = zeros(szx);
eU0   = zeros(szx); 
eV0   = zeros(szx);
OMEG  = zeros(szx); 
eOMEG = zeros(szx);

e = GSS.East;
n = GSS.North;
ve = GSS.EastVelocity;
vn = GSS.NorthVelocity;
eve = GSS.EastVelocityError;
evn = GSS.NorthVelocityError;
npe = GSS.IncludedPoints;

Iexcl = ~npe;
eve(Iexcl) = eve(Iexcl)*magnitf;
evn(Iexcl) = evn(Iexcl)*magnitf;

hwait = waitbar(0,'Strain calculation in progress, please wait...');
fprintf('\n\nStrain calculation in progress...\n');
for m = 1:szx(1)
    for mm = 1:szx(2)
        em = X(m,mm);
        nm = Y(m,mm);
        %------------------------------------------------------------------
        [Emax,Emin,Phi,eEmax,eEmin,n0,flc,fls,u0,v0,eu0,ev0,...
            omeg,eomeg] = gs_strain(em,nm,e,n,ve,vn,eve,evn,npe,d0,Options);
        %------------------------------------------------------------------
        EMAX(m,mm)  = Emax;  
        EMIN(m,mm)  = Emin; 
        PHI(m,mm)   = Phi;
        eEMAX(m,mm) = eEmax; 
        eEMIN(m,mm) = eEmin;
        numb(m,mm)  = n0;
        FLC(m,mm)   = flc;   
        FLS(m,mm)   = fls;
        U0(m,mm)    = u0;    
        V0(m,mm)    = v0;
        eU0(m,mm)   = eu0;   
        eV0(m,mm)   = ev0;
        OMEG(m,mm)  = omeg;  
        eOMEG(m,mm) = eomeg;
    end
    waitbar(m/szx(1),hwait);
end
delete(hwait);
fprintf('... Strain calculation completed\n\n');

GSSOut = GSS;
GSSOut.ScaleFactor = d0;
GSSOut.MaxStrainMatrix = EMAX;
GSSOut.MinStrainMatrix = EMIN;
GSSOut.AzimuthStrainMatrix = PHI;
GSSOut.MaxStrainMatrixError = eEMAX;
GSSOut.MinStrainMatrixError = eEMIN;
GSSOut.NumberMatrix = numb;
GSSOut.ComputationFlag = FLC;
GSSOut.SignificanceFlag = FLS;
GSSOut.UTranslation = U0;
GSSOut.VTranslation = V0;
GSSOut.UTranslationError = eU0;
GSSOut.VTranslationError = eV0;
GSSOut.Rotation = OMEG;
GSSOut.RotationError = eOMEG;